{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.gridspec import GridSpec\n",
    "plt.style.use('./deeplearning.mplstyle')\n",
    "import tensorflow as tf\n",
    "from tensorflow.keras.models import Sequential\n",
    "from tensorflow.keras.layers import Dense, LeakyReLU\n",
    "from tensorflow.keras.activations import linear, relu, sigmoid\n",
    "%matplotlib widget\n",
    "from matplotlib.widgets import Slider\n",
    "from lab_utils_common import dlc\n",
    "from autils import plt_act_trio\n",
    "from lab_utils_relu import *\n",
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=UserWarning)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"2\"></a>\n",
    "## 2 - ReLU Activation\n",
    "This week, a new activation was introduced, the Rectified Linear Unit (ReLU). \n",
    "$$ a = max(0,z) \\quad\\quad\\text {# ReLU function} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt_act_trio()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img align=\"right\" src=\"./images/C2_W2_ReLu.png\"     style=\" width:380px; padding: 10px 20px; \" >\n",
    "The example from the lecture on the right shows an application of the ReLU. In this example, the derived \"awareness\" feature is not binary but has a continuous range of values. The sigmoid is best for on/off or binary situations. The ReLU provides a continuous linear relationship. Additionally it has an 'off' range where the output is zero.     \n",
    "The \"off\" feature makes the ReLU a Non-Linear activation. Why is this needed? Let's examine this below. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Why Non-Linear Activations?  \n",
    "<img align=\"left\" src=\"./images/C2_W2_ReLU_Graph.png\"     style=\" width:250px; padding: 10px 20px; \" > The function shown is composed of linear pieces (piecewise linear). The slope is consistent during the linear portion and then changes abruptly at transition points. At transition points, a new linear function is added which, when added to the existing function, will produce the new slope. The new function is added at transition point but does not contribute to the output prior to that point. The Non-linear activation function is responsible for disabling the input prior to and sometimes after the transition points. The following exercise provides a more tangible example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The exercise will use the network below in a regression problem where you must model a piecewise linear target :\n",
    "<img align=\"center\" src=\"./images/C2_W2_ReLU_Network.png\"     style=\" width:650px; padding: 10px 20px; \">  \n",
    "The network has 3 units in the first layer. Each will be responsible for one piece of the target. Unit 0 is pre-programmed and fixed to map the the first segment. You will modify weights and biases in unit 1 and 2 to model the 2nd and 3rd segment. The output unit is also fixed and simply sums the outputs of the first layer.  \n",
    "\n",
    "Using the sliders below, modify weights and bias to match the target. \n",
    "Hints: Start with `w1` and `b1` and leave `w2` and `b2` zero until you match the 2nd segment. Clicking rather than sliding is quicker.  If you have trouble, don't worry, the text below will describe this in more detail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = plt_relu_ex()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " \n",
    "The goal of this exercise is to appreciate how the ReLU's non-linear behavior provides the needed ability to turn functions off until they are needed. Let's see how this worked in this example.\n",
    "<img align=\"right\" src=\"./images/C2_W2_ReLU_Plot.png\"     style=\" width:600px; padding: 10px 20px; \"> \n",
    "The plots on the right contain the output of the units in the first layer.   \n",
    "Starting at the top, unit 0 is responsible for the first segment marked with a 1. Both the linear function $z$ and the function following the ReLU $a$ are shown. You can see that the ReLU cuts off the function after the interval [0,1]. This is key as it prevents from interfering with the following segment. \n",
    "\n",
    "Unit 1 is responsible for the 2nd segment. Here the ReLU kept this unit quiet until after x is 1. Since the first unit is not contributing, the slope for unit 1, $w^{[1]}_1$, is just the slope of the target line. The bias must be adjusted to keep the output negative until x has reached 1.\n",
    "\n",
    "Unit 2 is responsible for the 3rd segment. The ReLU again zeros the output until x reaches the right value.The slope of the unit, $w^{[1]}_2$, must be set so that the sum of unit 1 and 2 have the desired slope. The bias is again adjusted to keep the output negative until x has reached 2.\n",
    "\n",
    "The \"off\" or disable feature  of the ReLU activation enables models to stitch together linear segments to model complex non-linear functions.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# =======End of new Material==="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## New Activation\n",
    "This week a new activation was introduced, the Rectified Linear Unit (ReLU). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plt_act_trio():\n",
    "    X = np.linspace(-5,5,100)\n",
    "    fig,ax = plt.subplots(1,3, figsize=(6,2))\n",
    "    widgvis(fig)\n",
    "    ax[0].plot(X,tf.keras.activations.linear(X))\n",
    "    ax[0].axvline(0, lw=0.3, c=\"black\")\n",
    "    ax[0].axhline(0, lw=0.3, c=\"black\")\n",
    "    ax[0].set_title(\"linear\")\n",
    "    ax[1].plot(X,tf.keras.activations.sigmoid(X))\n",
    "    ax[1].axvline(0, lw=0.3, c=\"black\")\n",
    "    ax[1].axhline(0, lw=0.3, c=\"black\")\n",
    "    ax[1].set_title(\"sigmoid\")\n",
    "    ax[2].plot(X,tf.keras.activations.relu(X))\n",
    "    ax[2].axhline(0, lw=0.3, c=\"black\")\n",
    "    ax[2].axvline(0, lw=0.3, c=\"black\")\n",
    "    ax[2].set_title(\"relu\")\n",
    "    fig.suptitle(\"Common Activation Functions\", fontsize=14)\n",
    "    fig.tight_layout(pad=0.2)\n",
    "    plt.show()\n",
    "plt_act_trio()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<right images>\n",
    "    The example on the right shows an application of the ReLu. In this example, the \"awareness\" feature is not binary but has a range that varies from 0 to a larger value. The sigmoid is best for on/off or binary situations. The ReLu provides a linear relationship And an 'off' range where the output is zero. The \"off\" feature makes the ReLu a Non-Linear activation. Why is this needed? Let's examine this with the example below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.linspace(0,2*np.pi, 100)\n",
    "y = np.cos(X)+1\n",
    "y[50:100]=0\n",
    "fig,ax = plt.subplots(1,1, figsize=(2,2))\n",
    "widgvis(fig)\n",
    "ax.plot(X,y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w10 = np.array([[-1]])\n",
    "b10 = np.array([2.6])\n",
    "d10 = Dense(1, activation = \"linear\", input_shape = (1,), weights=[w10,b10])\n",
    "z10 = d10(X.reshape(-1,1))\n",
    "a10 = relu(z10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plt_act1(y,z,a):\n",
    "    fig,ax = plt.subplots(1,3, figsize=(6,2.5))\n",
    "    widgvis(fig)\n",
    "    ax[0].plot(X,y,label=\"target\")\n",
    "    ax[0].axvline(0, lw=0.3, c=\"black\")\n",
    "    ax[0].axhline(0, lw=0.3, c=\"black\")\n",
    "    ax[0].set_title(\"y - target\")\n",
    "    ax[1].plot(X,y, label=\"target\")\n",
    "    ax[1].plot(X,z, c=dlc[\"dldarkred\"],label=\"z\")\n",
    "    ax[1].axvline(0, lw=0.3, c=\"black\")\n",
    "    ax[1].axhline(0, lw=0.3, c=\"black\")\n",
    "    ax[1].set_title(\"z = wX+b\")\n",
    "    ax[1].legend(loc=\"upper center\")\n",
    "    ax[2].plot(X,y, label=\"target\")\n",
    "    ax[2].plot(X,a, c=dlc[\"dldarkred\"],label=\"ReLu(z)\")\n",
    "    ax[2].axhline(0, lw=0.3, c=\"black\")\n",
    "    ax[2].axvline(0, lw=0.3, c=\"black\")\n",
    "    ax[2].set_title(\"with relu\")\n",
    "    ax[2].legend()\n",
    "    fig.suptitle(\"Role of Activation\", fontsize=14)\n",
    "    fig.tight_layout(pad=0.2)\n",
    "    return(ax)\n",
    "\n",
    "\n",
    "def plt_add_notation(ax):\n",
    "    ax[1].annotate(text = \"matches\\n here\", xy =(1.5,1.0), \n",
    "                   xytext = (0.1,-1.5), fontsize=10,\n",
    "                  arrowprops=dict(facecolor=dlc[\"dlpurple\"],width=2, headwidth=8))\n",
    "    ax[1].annotate(text = \"but not\\n here\", xy =(5,-2.5), \n",
    "                   xytext = (1,-3), fontsize=10,\n",
    "                  arrowprops=dict(facecolor=dlc[\"dlpurple\"],width=2, headwidth=8))\n",
    "    ax[2].annotate(text = \"ReLu\\n 'off'\", xy =(2.6,0), \n",
    "                   xytext = (0.1,0.1), fontsize=10,\n",
    "                  arrowprops=dict(facecolor=dlc[\"dlpurple\"],width=2, headwidth=8))\n",
    "    \n",
    "ax = plt_act1(y,z10,a10)\n",
    "plt_add_notation(ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.linspace(0,2*np.pi, 100)\n",
    "y = np.cos(X)+1\n",
    "y[0:49]=0\n",
    "fig,ax = plt.subplots(1,1, figsize=(2,2))\n",
    "widgvis(fig)\n",
    "ax.plot(X,y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w11 = np.array([[1]])\n",
    "b11 = np.array([-3.7])\n",
    "d11 = Dense(1, activation = \"linear\", input_shape = (1,), weights=[w11,b11])\n",
    "z11 = d11(X.reshape(-1,1))\n",
    "a11 = relu(z11)\n",
    "plt_act1(y,z11,a11)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.linspace(0,2*np.pi, 100)\n",
    "y = np.cos(X)+1\n",
    "X=X.reshape(-1,1)\n",
    "yhat = relu(d10(X)) + relu(d11(X))\n",
    "\n",
    "fig,ax = plt.subplots(1,2, figsize=(4,2))\n",
    "widgvis(fig)\n",
    "ax[0].plot(X,y)\n",
    "ax[1].plot(X,y)\n",
    "ax[1].plot(X,yhat)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X=X.reshape(-1,1)\n",
    "yhat = relu(d10(X)) + relu(d11(X))\n",
    "fig,ax = plt.subplots(1,1, figsize=(2,2))\n",
    "widgvis(fig)\n",
    "ax.plot(X,y)\n",
    "ax.plot(X,yhat)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yhat = model.predict(X.reshape(-1,1))\n",
    "fig,ax = plt.subplots(1,1, figsize=(2,2))\n",
    "widgvis(fig)\n",
    "ax.plot(X,y)\n",
    "ax.plot(X,yhat)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "model = Sequential(\n",
    "    [ \n",
    "        d10,\n",
    "        tf.layers.Activation(activations.relu),\n",
    "        d11,\n",
    "        tf.layers.Activation(activations.relu),\n",
    "        Dense(1, activation = 'linear')\n",
    "    ]\n",
    ")\n",
    "model.compile(\n",
    "    loss=tf.keras.losses.MeanSquaredError(),\n",
    "    optimizer=tf.keras.optimizers.Adam(0.1),\n",
    ")\n",
    "\n",
    "model.fit(\n",
    "    X,y,\n",
    "    epochs=1000\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "model = Sequential(\n",
    "    [ \n",
    "        Dense(1,activation=\"relu\", name = 'l1'),\n",
    "        Dense(1,activation=\"linear\", name = 'l2')\n",
    "    ]\n",
    ")\n",
    "model.compile(\n",
    "    loss=tf.keras.losses.MeanSquaredError(),\n",
    "    optimizer=tf.keras.optimizers.Adam(0.01),\n",
    ")\n",
    "\n",
    "model.fit(\n",
    "    X,y,\n",
    "    epochs=10\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yhat = model.predict(X)\n",
    "yhat[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(1,1, figsize=(4,4))\n",
    "ax.plot(X,y)\n",
    "ax.plot(X,yhat, c=dlc[\"dldarkred\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l1 = model.get_layer('l1')\n",
    "l2 = model.get_layer('l2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l1.get_weights()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l2.get_weights()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l1 = model.get_layer('l1')\n",
    "l2 = model.get_layer('l2')\n",
    "l1.get_weights()\n",
    "l2.get_weights()\n",
    "\n",
    "w1 = np.array([[-1]])\n",
    "b1 = np.array([1])\n",
    "l1.set_weights([w1,b1])\n",
    "\n",
    "w2 = np.array([[1]])\n",
    "b2 = np.array([0])\n",
    "l2.set_weights([w2,b2])\n",
    "model.fit(\n",
    "    X,y,\n",
    "    epochs=100\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l2.set_weights([w2,b2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yhat = model.predict(X)\n",
    "fig,ax = plt.subplots(1,1, figsize=(4,4))\n",
    "ax.plot(X,y)\n",
    "ax.plot(X,yhat, c=dlc[\"dldarkred\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2-D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import warnings\n",
    "from matplotlib import cm\n",
    "from matplotlib.patches import FancyArrowPatch\n",
    "from matplotlib.colors import ListedColormap, LinearSegmentedColormap\n",
    "import matplotlib.colors as colors\n",
    "from lab_utils_common import dlc\n",
    "dkcolors = plt.cm.Paired((1,3,7,9,5,11))\n",
    "ltcolors = plt.cm.Paired((0,2,6,8,4,10))\n",
    "dkcolors_map = mpl.colors.ListedColormap(dkcolors)\n",
    "ltcolors_map = mpl.colors.ListedColormap(ltcolors)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plt_mc_data(ax, X, y, classes,  class_labels=None, map=plt.cm.Paired, \n",
    "                legend=False, size=50, m='o', equal_xy = False):\n",
    "    \"\"\" Plot multiclass data. Note, if equal_xy is True, setting ylim on the plot may not work \"\"\"\n",
    "    for i in range(classes):\n",
    "        idx = np.where(y == i)\n",
    "        col = len(idx[0])*[i]\n",
    "        label = class_labels[i] if class_labels else \"c{}\".format(i)\n",
    "        ax.scatter(X[idx, 0], X[idx, 1],  marker=m,\n",
    "                    c=col, vmin=0, vmax=map.N, cmap=map,\n",
    "                    s=size, label=label)\n",
    "    if legend: ax.legend()\n",
    "    if equal_xy: ax.axis(\"equal\")\n",
    "\n",
    "def plt_mc(X_train,y_train,classes):\n",
    "    css = np.unique(y_train)\n",
    "    fig,ax = plt.subplots(1,1,figsize=(3,3))\n",
    "    fig.canvas.toolbar_visible = False\n",
    "    fig.canvas.header_visible = False\n",
    "    fig.canvas.footer_visible = False\n",
    "    plt_mc_data(ax, X_train,y_train,classes, map=dkcolors_map, legend=True, size=10, equal_xy = False)\n",
    "    ax.set_title(\"Multiclass Data\")\n",
    "    ax.set_xlabel(\"x0\")\n",
    "    ax.set_ylabel(\"x1\")\n",
    "    return(ax)\n",
    "    \n",
    "def plot_cat_decision_boundary_mc(ax, X, predict , class_labels=None, legend=False, vector=True):\n",
    "\n",
    "    # create a mesh to points to plot\n",
    "    x_min, x_max = X[:, 0].min(), X[:, 0].max()\n",
    "    y_min, y_max = X[:, 1].min(), X[:, 1].max()\n",
    "    h = max(x_max-x_min, y_max-y_min)/200\n",
    "    xx, yy = np.meshgrid(np.arange(x_min, x_max, h),\n",
    "                         np.arange(y_min, y_max, h))\n",
    "    points = np.c_[xx.ravel(), yy.ravel()]\n",
    "    #print(\"points\", points.shape)\n",
    "    #print(\"xx.shape\", xx.shape)\n",
    "\n",
    "    #make predictions for each point in mesh\n",
    "    if vector:\n",
    "        Z = predict(points)\n",
    "    else:\n",
    "        Z = np.zeros((len(points),))\n",
    "        for i in range(len(points)):\n",
    "            Z[i] = predict(points[i].reshape(1,2))\n",
    "    Z = Z.reshape(xx.shape)\n",
    "\n",
    "    #contour plot highlights boundaries between values - classes in this case\n",
    "    ax.contour(xx, yy, Z, linewidths=1) \n",
    "    #ax.axis('tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.random.rand(300, 2)\n",
    "y = np.sqrt( X[:,0]**2 + X[:,1]**2 ) < 0.6\n",
    "#y = np.logical_and( X[:,0] < 0.5, X[:,1] < 0.5 ).astype(int)\n",
    "y.shape\n",
    "plt_mc(X,y,2,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "model = Sequential(\n",
    "    [ \n",
    "        Dense(2,activation=\"relu\",    name = 'l1'),\n",
    "        Dense(1,activation=\"sigmoid\", name = 'l2')\n",
    "    ]\n",
    ")\n",
    "\n",
    "model.compile(\n",
    "    loss=tf.keras.losses.MeanSquaredError(),\n",
    "    optimizer=tf.keras.optimizers.Adam(0.01),\n",
    ")\n",
    "\n",
    "model.fit(\n",
    "    X,y,\n",
    "    epochs=150\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt_mc(X,y,2,)\n",
    "predict = lambda  x: (model.predict(x) > 0.5).astype(int)\n",
    "plot_cat_decision_boundary_mc(ax, X, predict, legend = True, vector=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l1 = model.get_layer(\"l1\")\n",
    "W1,b1 = l1.get_weights()\n",
    "l2 = model.get_layer(\"l2\")\n",
    "W2,b2 = l2.get_weights()\n",
    "print(W1,b1)\n",
    "print(W2,b2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = np.array([0.4,0.60])\n",
    "np.dot( np.dot(x0,W1) + b1, W2) + b2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn import cluster, datasets, mixture\n",
    "from sklearn.neighbors import kneighbors_graph\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from itertools import cycle, islice\n",
    "\n",
    "np.random.seed(0)\n",
    "\n",
    "# ============\n",
    "# Generate datasets. We choose the size big enough to see the scalability\n",
    "# of the algorithms, but not too big to avoid too long running times\n",
    "# ============\n",
    "n_samples = 500\n",
    "noisy_circles = datasets.make_circles(n_samples=n_samples, factor=0.5, noise=0.05)\n",
    "noisy_moons = datasets.make_moons(n_samples=n_samples, noise=0.05)\n",
    "blobs = datasets.make_blobs(n_samples=n_samples, random_state=8)\n",
    "no_structure = np.random.rand(n_samples, 2), None\n",
    "\n",
    "# Anisotropicly distributed data\n",
    "random_state = 170\n",
    "X, y = datasets.make_blobs(n_samples=n_samples, random_state=random_state)\n",
    "transformation = [[0.6, -0.6], [-0.4, 0.8]]\n",
    "X_aniso = np.dot(X, transformation)\n",
    "aniso = (X_aniso, y)\n",
    "\n",
    "# blobs with varied variances\n",
    "varied = datasets.make_blobs(\n",
    "    n_samples=n_samples, cluster_std=[1.0, 2.5, 0.5], random_state=random_state\n",
    ")\n",
    "\n",
    "# ============\n",
    "# Set up cluster parameters\n",
    "# ============\n",
    "plt.figure(figsize=(9 * 2 + 3, 13))\n",
    "plt.subplots_adjust(\n",
    "    left=0.02, right=0.98, bottom=0.001, top=0.95, wspace=0.05, hspace=0.01\n",
    ")\n",
    "\n",
    "plot_num = 1\n",
    "\n",
    "default_base = {\n",
    "    \"quantile\": 0.3,\n",
    "    \"eps\": 0.3,\n",
    "    \"damping\": 0.9,\n",
    "    \"preference\": -200,\n",
    "    \"n_neighbors\": 3,\n",
    "    \"n_clusters\": 3,\n",
    "    \"min_samples\": 7,\n",
    "    \"xi\": 0.05,\n",
    "    \"min_cluster_size\": 0.1,\n",
    "}\n",
    "\n",
    "datasets = [\n",
    "    (\n",
    "        noisy_circles,\n",
    "        {\n",
    "            \"damping\": 0.77,\n",
    "            \"preference\": -240,\n",
    "            \"quantile\": 0.2,\n",
    "            \"n_clusters\": 2,\n",
    "            \"min_samples\": 7,\n",
    "            \"xi\": 0.08,\n",
    "        },\n",
    "    ),\n",
    "    (\n",
    "        noisy_moons,\n",
    "        {\n",
    "            \"damping\": 0.75,\n",
    "            \"preference\": -220,\n",
    "            \"n_clusters\": 2,\n",
    "            \"min_samples\": 7,\n",
    "            \"xi\": 0.1,\n",
    "        },\n",
    "    ),\n",
    "    (\n",
    "        varied,\n",
    "        {\n",
    "            \"eps\": 0.18,\n",
    "            \"n_neighbors\": 2,\n",
    "            \"min_samples\": 7,\n",
    "            \"xi\": 0.01,\n",
    "            \"min_cluster_size\": 0.2,\n",
    "        },\n",
    "    ),\n",
    "    (\n",
    "        aniso,\n",
    "        {\n",
    "            \"eps\": 0.15,\n",
    "            \"n_neighbors\": 2,\n",
    "            \"min_samples\": 7,\n",
    "            \"xi\": 0.1,\n",
    "            \"min_cluster_size\": 0.2,\n",
    "        },\n",
    "    ),\n",
    "    (blobs, {\"min_samples\": 7, \"xi\": 0.1, \"min_cluster_size\": 0.2}),\n",
    "    (no_structure, {}),\n",
    "]\n",
    "\n",
    "datasets = [\n",
    "    (no_structure, {}),\n",
    "]\n",
    "\n",
    "for i_dataset, (dataset, algo_params) in enumerate(datasets):\n",
    "    # update parameters with dataset-specific values\n",
    "    params = default_base.copy()\n",
    "    params.update(algo_params)\n",
    "\n",
    "    X, y = dataset\n",
    "\n",
    "    # normalize dataset for easier parameter selection\n",
    "    X = StandardScaler().fit_transform(X)\n",
    "\n",
    "    # estimate bandwidth for mean shift\n",
    "    bandwidth = cluster.estimate_bandwidth(X, quantile=params[\"quantile\"])\n",
    "\n",
    "    # connectivity matrix for structured Ward\n",
    "    connectivity = kneighbors_graph(\n",
    "        X, n_neighbors=params[\"n_neighbors\"], include_self=False\n",
    "    )\n",
    "    # make connectivity symmetric\n",
    "    connectivity = 0.5 * (connectivity + connectivity.T)\n",
    "\n",
    "    # ============\n",
    "    # Create cluster objects\n",
    "    # ============\n",
    "    ms = cluster.MeanShift(bandwidth=bandwidth, bin_seeding=True)\n",
    "    two_means = cluster.MiniBatchKMeans(n_clusters=params[\"n_clusters\"])\n",
    "    ward = cluster.AgglomerativeClustering(\n",
    "        n_clusters=params[\"n_clusters\"], linkage=\"ward\", connectivity=connectivity\n",
    "    )\n",
    "    spectral = cluster.SpectralClustering(\n",
    "        n_clusters=params[\"n_clusters\"],\n",
    "        eigen_solver=\"arpack\",\n",
    "        affinity=\"nearest_neighbors\",\n",
    "    )\n",
    "    dbscan = cluster.DBSCAN(eps=params[\"eps\"])\n",
    "    optics = cluster.OPTICS(\n",
    "        min_samples=params[\"min_samples\"],\n",
    "        xi=params[\"xi\"],\n",
    "        min_cluster_size=params[\"min_cluster_size\"],\n",
    "    )\n",
    "    affinity_propagation = cluster.AffinityPropagation(\n",
    "        damping=params[\"damping\"], preference=params[\"preference\"], random_state=0\n",
    "    )\n",
    "    average_linkage = cluster.AgglomerativeClustering(\n",
    "        linkage=\"average\",\n",
    "        affinity=\"cityblock\",\n",
    "        n_clusters=params[\"n_clusters\"],\n",
    "        connectivity=connectivity,\n",
    "    )\n",
    "    birch = cluster.Birch(n_clusters=params[\"n_clusters\"])\n",
    "    gmm = mixture.GaussianMixture(\n",
    "        n_components=params[\"n_clusters\"], covariance_type=\"full\"\n",
    "    )\n",
    "\n",
    "    clustering_algorithms = (\n",
    "        (\"MiniBatch\\nKMeans\", two_means),\n",
    "        (\"Affinity\\nPropagation\", affinity_propagation),\n",
    "        (\"MeanShift\", ms),\n",
    "        (\"Spectral\\nClustering\", spectral),\n",
    "        (\"Ward\", ward),\n",
    "        (\"Agglomerative\\nClustering\", average_linkage),\n",
    "        (\"DBSCAN\", dbscan),\n",
    "        (\"OPTICS\", optics),\n",
    "        (\"BIRCH\", birch),\n",
    "        (\"Gaussian\\nMixture\", gmm),\n",
    "    )\n",
    "\n",
    "    for name, algorithm in clustering_algorithms:\n",
    "        t0 = time.time()\n",
    "\n",
    "        # catch warnings related to kneighbors_graph\n",
    "        with warnings.catch_warnings():\n",
    "            warnings.filterwarnings(\n",
    "                \"ignore\",\n",
    "                message=\"the number of connected components of the \"\n",
    "                + \"connectivity matrix is [0-9]{1,2}\"\n",
    "                + \" > 1. Completing it to avoid stopping the tree early.\",\n",
    "                category=UserWarning,\n",
    "            )\n",
    "            warnings.filterwarnings(\n",
    "                \"ignore\",\n",
    "                message=\"Graph is not fully connected, spectral embedding\"\n",
    "                + \" may not work as expected.\",\n",
    "                category=UserWarning,\n",
    "            )\n",
    "            print(X.shape,algorithm)\n",
    "            algorithm.fit(X)\n",
    "\n",
    "        t1 = time.time()\n",
    "        if hasattr(algorithm, \"labels_\"):\n",
    "            y_pred = algorithm.labels_.astype(int)\n",
    "        else:\n",
    "            y_pred = algorithm.predict(X)\n",
    "\n",
    "        plt.subplot(len(datasets), len(clustering_algorithms), plot_num)\n",
    "        if i_dataset == 0:\n",
    "            plt.title(name, size=18)\n",
    "\n",
    "        colors = np.array(\n",
    "            list(\n",
    "                islice(\n",
    "                    cycle(\n",
    "                        [\n",
    "                            \"#377eb8\",\n",
    "                            \"#ff7f00\",\n",
    "                            \"#4daf4a\",\n",
    "                            \"#f781bf\",\n",
    "                            \"#a65628\",\n",
    "                            \"#984ea3\",\n",
    "                            \"#999999\",\n",
    "                            \"#e41a1c\",\n",
    "                            \"#dede00\",\n",
    "                        ]\n",
    "                    ),\n",
    "                    int(max(y_pred) + 1),\n",
    "                )\n",
    "            )\n",
    "        )\n",
    "        # add black color for outliers (if any)\n",
    "        colors = np.append(colors, [\"#000000\"])\n",
    "        plt.scatter(X[:, 0], X[:, 1], s=10, color=colors[y_pred])\n",
    "\n",
    "        plt.xlim(-2.5, 2.5)\n",
    "        plt.ylim(-2.5, 2.5)\n",
    "        plt.xticks(())\n",
    "        plt.yticks(())\n",
    "        plt.text(\n",
    "            0.99,\n",
    "            0.01,\n",
    "            (\"%.2fs\" % (t1 - t0)).lstrip(\"0\"),\n",
    "            transform=plt.gca().transAxes,\n",
    "            size=15,\n",
    "            horizontalalignment=\"right\",\n",
    "        )\n",
    "        plot_num += 1\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
